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A new phenomenological model of turbulent fluctuations is constructed by considering the La- 
grangian dynamics of 4 points (the tetrad). The closure of the equations of motion is achieved by 
^T *' postulating an anisotropic, i.e. tetrad shape dependent, relation of the local pressure and the velocity 

gradient defined on the tetrad. The non-local contribution to the pressure and the incoherent small 
scale fluctuations are modeled as Gaussian white "noise". The resulting stochastic model for the 
coarse-grained velocity gradient is analyzed approximately, yielding predictions for the probability 
distribution functions of different 2nd and 3d order invariants. The results are compared with the 
direct numerical simulation of the Navier-Stokes. The model provides a reasonable representation 
of the non-linear dynamics involved in energy transfer and vortex stretching and allows to study 
^ . interesting aspects of the statistical geometry of turbulence, e.g. vorticity /strain alignment. In a 

' state with a constant energy flux (and K41 power spectrum) it exhibits the anomalous scaling of 

fS) , high moments associated with formation of high gradient sheets - events associated with large energy 

■ transfer. An approach to the more complete analysis of the stochastic model, properly including 

' the effect of fluctuations, is outlined and will enable further quantitative juxtaposition of the model 

O I with the results of the DNS. 

a^ 

I. INTRODUCTION 

I ' The old problem of hydrodynamic turbulence has in recent years attracted resurgent interest stimulated by the 
O I new generation of laboratory experiments and the newly acquired ability of the direct numerical simulations to 
probe interesting aspects of turbulence. In the light of the new ideas and developments there has also been new 
appreciation of the seminal contributions of Kolmogorov, reviewed in a recent book by Frisch jl], and of Kraichnan 
to whom the present volume and this article are dedicated. The key issues and the progress of the last years have 
' been well reviewed and are well represented in the present Festschrift. Much effort has been dedicated to 1) 
,: documenting and understanding the anomalous (i.e. non-Kolmogorov 41) scaling of high moments associated 
^ with intermittency, 2) understanding the structure and the local geometry of the intermittent regions of the flow [j7[-p5[ . 
5^ ^ On the theory side new ideas derived from the new understanding of anomalous scaling of the Passive Scalar problem 
p6| pO| and of the Burger's turbulence P,pT|-p3[ both pioneered by Kraichnan |2^]. Yet, the theoretical description 
of turbulence based on first principles, i.e. on a controlled approximation to the Navier-Stokes equations, is still over 
the horizon and to proceed in the right direction one must rely on phenomenology. One reason for pursuing the 
modeling approach is the need to bridge the existing gulf between our understanding of the scaling of turbulent 
fluctuations and their structure or "statistical geometry" [^^. A step in this direction will be the subject of the 
present paper. 

Our goal here is to advance a phenomenological model for the probability distribution function (PDF) of turbulent 
velocity fluctuations. We shall start by noting that the longitudinal velocity difference between two observation 
points p6|| , while being most readily observable, seems a poor candidate for a fundamental dynamical field in terms 
of which to attempt a closed statistical description. The intuitive reason is that the longitudinal velocity difference 
senses only one of the 8 locally independent components of the velocity gradient tensor which govern the dynamics 
of the velocity field. Instead we shall choose the fundamental field to be the coarse grained velocity gradient tensor 
Mab = /p drdaVb{r) defined over a region T with characteristic scale R lying in the inertial range. This region may be 
best thought of as a local correlation volume of the velocity gradient coarse-grained on scale i? - an "eddy" of sort. 
The phenomenological model then will be based on the Lagrangian dynamics of the F-volume, parametrized by four 
points - the "tetrad" - and its strain and vorticity flelds as described by Mab- Effort will be made to preserve the 
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essential non-linear dynamics governing evolution of coarse grained strain and vorticity and the concomitant distortion 
of the Lagrangian volume. This dynamics expresses the fundamental constraints due to the conservation of energy 
and circulation. In contrast, the dynamics of velocity fluctuations arising from the scales smaller than that of the 
tetrad (and generating incoherent motion of the points) will be modeled as a Gaussian white process obeying K41 
scaling. Essential element of the theory advanced below will be the decomposition of the pressure into the local part 
determined by the M-field via incompressibility and the non-local part due to the contribution of distant regions, 
which will again be modeled as a Gaussian random force. Such an approximation for the pressure may perhaps be 
justifiable in large spatial dimension [^7|- ^ , p^ , but we shall be content with exploring its consequences and comparing 
the results with the direct numerical simulations. 

We cannot hope to review here the evolution of the phenomenological modeling ideas, yet we shall put the present 
work into the context of two recent efforts. The "PDF models" of Pope and coworkers attempt to close the 
equation for the velocity probability distribution function (PDF) on the level of one point: in contrast, our model 
deals with relative velocity differences on 4 points which naturally brings in Kolmogorov's ideas and allows to address 
the intermittency phenomenon. The two approaches however share the need to model the pressure Hessian/strain- 
rate correlations (in our case on the coarse-grained level) and share the realization that this model is improved by 
incorporating dynamical information about local anisotropy (in our case furnished by the moment of inertia 

tensor of the evolving Lagrangian volume) . Another point of reference is the two — point PDF closure advanced by 
Yakhot 1 29 1 on the basis of the work of Kraichnan and Polyakov |2^ . There too one arrived at a Fokker-Planck 
type equation for the PDF of velocity differences at given point separation, yet the approach differs from the present 
one in the treatment of the correlations of large and small scale fluctuations and our approach, by virtue of tracking 
a tetrad rather than a pair will retain more of the geometry of the flow. 

The model will be presented in Section 2 in the form of the stochastic equations of motion for two tensors specifying 
the coarse-grained velocity gradient and the shape (i.e. moment of inertia) of the evolving Lagrangian volume. We 
shall write down the corresponding Fokker-Planck equation for the Probability Distribution Function and discuss 
the energy transfer considerations which played the key role in the formulation of the model. Section 3 relates the 
deterministic aspects of the model to the so called Restricted Euler (RE) dynamics that has been investigated by 
Vieillefosse originally in the context of the finite time singularities (see also Leorat Q|) and Cantwell and coworkers 
with the emphasis on the local topology of the flow. RE describes the evolution of the velocity gradient at a 
point within an isotropic approximation for the pressure which allows to close the Euler equation locally. We shall 
see that elevation of the dynamics to the coarse-grained level and the introduction of the second dynamical field to 
keep track of the shape of the Lagrangian correlation volume (which depends on the history of the strain) allows 
to go beyond the isotropic pressure approximation: the unphysical finite time singularity of RE is removed, while 
the sensible short time dynamical properties (related for example to the vorticity strain alignment Q) are retained. 
Finally in Section 3 the deterministic dynamics will be compared with the empirical " mean field" equation of motion 
for the coarse grained velocity tensor M, constructed from the conditional average {^M|M^ measured in the DNS of 

the Navier-Stokes at R\ = 85. In Section 4 we return to the stochastic model and write down the formal solution of 
the Fokker-Planck equation in terms of the path integral relating the probability of a given coarse grained velocity 
gradient on a given inertial range tetrad to the velocity PDF on the integral scale. This path integral representation 
serves as a point of departure for the semiclassical approximation. It also has a well defined deterministic limit where 
the effect of the stochastic terms in the tetrad dynamics can be neglected. In the latter limit the probability of 
"observing" any given coarse-grained velocity gradient on an inertial range tetrad is determined by the probability 
of its integral scale preimage and can be calculated by integrating the equations of motion backward in time. This 
crude but simple approximation is employed in Section 5 in order to gain insight into the behavior of the model and 
to elaborate its statistical predictions emphasizing energy transfer, enstrophy and its production and the alignment 
of vorticity and strain. The comparison of the results with the direct numerical simulation of the Navier-Stokes 
equations is quite encouraging. The calculated probability distribution function also exhibits anomalous scaling of 
high moments. In conclusion, Section 6 is a summary and the outline of further inquiry. 



II. THE MINIMAL MODEL 



The minimal parametrization of the T volume is a tetrahedron (more generally a d-dimensional simplex) defined 
by four, hence tetrad, (or d+1) Lagrangian points, ra{t), which upon elimination of the center of mass define a triad 
of vectors p^: pi ~ (ri — r2)/V2, p2 = {ri + r 2 — 2r3)/V6, = (fi + r 2 + — 3r4)/\/T2. It will be useful to treat 
this triad of vectors as 3 x 3 matrix, p° , where a is the spatial index. Analogously, by eliminating the center of mass 
velocity from the instantaneous velocity of the vertices, r^, one can define a triad of relative velocities, Vi. The coarse 
grained gradient field can now be defined simply by interpolation: 
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Mat,^{p-^)fv1-^-ftr{p-\). (1) 

Alternatively and more generally one may decompose the "observed" velocity differences, vf into a slow component 
arising from the scales greater than the radius of gyration, >R, represented by the coarse grained velocity gradient 
matrix Mab and the rapidly fluctuating incoherent component, Uj arising from scales ^ R. 

vt = p\Mab + ul (2) 

The strategy will be to derive the dynamics of Mab and while treating u as a Gaussian white noise with the 
statistics depending on the Kolmogorov's energy dissipation rate e as well as the instantaneous p and M. Lagrangian 
dynamics is governed by -^v^ = — Vpi + fi and -§iPi = v^, where is the random external force and Vpi is the 
pressure gradient, both properly differenced between the observation points. The dynamics of the coarse grained 
velocity gradient tensor, M, and the tetrad tensor p has the form: 



d 
di 



Mab + Mafc - ^abtrM^ = U, (3a) 

- p^Mba = <, (3b) 
Uab = Kk'^/trkk\ (3c) 



with matrix k = p~^. The left hand sides of (|3a,3b) describe the self-advection and stretching of the tetrad by the 
coherent (on the scale of p ) component of the velocity field. The right hand side of ( pa[ ) derives from the pressure 
gradient and the random force differences as well as from the coupling to incoherent small scale fluctuations. The 
Hab trM"^ term, with trYl = 1 on the left hand side represents the " local" component of the pressure needed to insure 
conservation of trM = as required by incompressibility. Tensor Hab is a measure of tetrad anisotropy representing 
the anisotropy of the Lagrangian correlation volume built up by prior evolution . This choice of the local term, 
in contrast with the simpler, isotropic form psf is dictated by the requirement that the pressure forces should do no 
work and drop out of the energy balance (see below). In addition it reintroduces proper reduction of the deterministic 
dynamics (left hand side of (Ba )) to 2D: i.e. a 2D velocity gradient configuration remains 2D provided that the tetrad 
has the shape of a filament, i.e. p^p is a tensor of rank one. The remaining non-local part of the pressure is subsumed 
in £,ab- 

We now define the stochastic components of the model appearing on the right hand sides of (3a-b). It is appealing 
to model the non-local part of the pressure retained in C, along with the contribution of small scales, as 5-correlated 
Gaussian random noise with the variance depending on the energy flux e as well as the instantaneous p and M. Let 
us consider a polynomial : 

i = ri + C,m + a{M^ -litrM^). (4) 

where 77 is a random matrix, and C a random function. The dimension of both ^ and is time^^ so that a is a 
constant. The last term is clearly not the most general one can write. It is however the one suggested by the numerical 



study of the statistics of the right hand side of (paj), originally by Borue and Orszag |39|. According to the DNS the 
average conditioned on M is not zero, but is reasonably well approximated by a(M^ — IltrM^) with < a < .8, 
depending on the scale, at least for the isotropic tetrad (i.e. 11 — 1/3). We shall assume that in the inertial range 
a is constant (which corresponds to keeping only the deterministic component of the 3d term in (4)) and take it to 
be a model parameter. This a term " renormalizes" the time scale of the deterministic dynamics described by the left 
hand side of the equation (|3a| ) and will have an important consequence for the energy transfer in the model, as we 
shall see shortly. 
Let 

2C e 1 

{Vab{t)-ncd{^)) = -^^7^ [Sackd - -^6ab5cd] ^{1) (5a) 

which is the simplest form respecting incompressibility, with Kolmogorov's energy flux e and C,, - a dimensionless 
parameter. Random r/ causes diffusion in velocity space; note that e has the dimension of corresponding diffusivity. 
The appearence of e in (5a) is further supported by the fact that rj fluctuations contribute to the energy transfer as 
we shall see below. In the " minimal" model, which we are now constructing we shall drop the possible multiplicative 
random field ( (see Eq.(^)). 

The fiuctuations of the small scale |M u can be resolved into parts which are longitudinal and transverse to pf: 
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«(<)w5(0)) = 2q|VirMMt + 2CWtrMM^ [p'S'^^S^, - p^p)) 5{t) (5b) 



where with the K41 scaling in mind we take the characteristic time to be the "eddy turnover" time, l/virMMJ. 
The longitudinal part of u in the p-equation (3b) would by itself produce Richardson diffusion behavior, (p^{t)) ^ et^ 



provided that the Kolmogorov scaling V irMMt ^e^/^p-^/a holds. However, Richardson diffusion would also arise 
from the non-Gaussian coherent stretching term pM, and the Gaussian longitudinal fluctuations, C||, do not appear 
to be essential. We shall set C|| =0. The transverse fluctuations C± however are very important, because in their 
absence the effect of volume preserving coherent stretching would lead to the rapid growth of anisotropy of the tetrad. 
The incoherent transverse velocity fluctuations act to redistribute the vertices of the tetrad uniformly on the surface 
of the p^ = trp^ p hyper-sphere in = 9 dimensions, thus introducing the isotropisation mechanism. The competition 
of the coherent stretching which leads to the growth of the radius of gyration p^ (both forward and backward in time) 
and the isotropisation over the p^ shell will play the key role in setting up the energy flux. 

The stochastic tetrad dynamics defined by Eq. (3-5) determines the Lagrangian transition probability from tetrad 
(M', p') to (M, p) at a time t later, G't(M, p|M', p'), which satisfies a Fokker-Planck equation 



Gt(M,p|M',p') = 5{M-M')5{p- p') 



- L ) Gt(M, p|M', p') = 5{M - M') 5{p - p') (6) 



with the evolution operator 



L = (1 - (iC - n.i trM?) - A plM.t (7) 

The invariant joint distribution, P(M, p), satisfying 

atP(M,p) = LP(M,p) = (8) 

can be interpreted as the Eulerian PDF P(M, p) provided that the normalization / dM P(M, p) = f is imposed [ p2[ . 
Equations (7-8) once supplemented with the boundary condition specifying the Eulerian PDF on the integral scale, 
p^ — L^, completely define our model. 

Before proceeding with the analysis of the Fokker-Planck equation (8), let us examine the energy balance, which 
was one of the key consideration in the formulation of the model: 



^dt < trVVt >p = J dM tr(pMMtp^) dtPiM, p) 



= -i^<V: tr{VV^) >p+a< tr{VV^M) >p 

+ ^Cr,e-C^Da + CU, (9) 

which is obtained by multiplying eq (||) by trVV^ (where = p\Mf,a) and averaging with respect to M. Note that 
the average < ... >p is taken at fixed p and remains a function of it. The first two terms on the right hand side may 
be identified as the divergence of the large scale energy flux and the eddy damping respectively. 

Note that the term originating from the deterministic component of the pressure Ilafc trM^ drops out: the particular 
form of flab was chosen for that purpose on the grounds that the pressure gradients should not contribute to the energy 
transfer as seen in the von Karman-Howarth derivation p3| . However, since is only the coarse grained and not the 
full local velocity, in contrast with the von Karman - Howarth analysis the divergence of the energy flux is balanced 
not directly by the viscous dissipation term, but by the eddy damping. There are also additional contributions due 
to the coupling with small scale fluctuations represented by the last three terms in (|^). The C,, term represents the 
diffusive component of the energy flux arising from the small scale fluctuations and the coupling of the tetrad to the 
neighboring regions (entering via Gaussian 77). The C_l/ term represents the transverse energy flux with: 



/VtrMMt irVvA - 2p5 /VtrMMt (MM.'^Vc 
op'] \ / p ■' \ 



(10) 
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which redistributes the energy within the = const shell, while the C±Dd is the diffusive contribution to the eddy 
damping: 



Dd = 18 ^VtrMMt tr [( pV " 3 P^) MM^^]^ p. (11a) 
This _D(i is reminiscent of the Smagorinsky form of eddy damping (popular in sub-grid simu lations [^ , ^ ) 

h 



but with significant difference that in the latter the c?-wave projector (pV ^ ^P^) appearing in (11a) is replaced by 



the simple scale factor p^. Hence, in contrast with the Smagorinsky model, our diffusive damping term is only active 
to the extent that the (vrrMMt (MM^)af, \ p tensor is correlated with {p^ p)ab tensor within the p^ shell. Strictly 



speaking Dd is not positive definite and its interpretation as the damping term is contingent on the expectation that 
the tetrad dynamics builds up the alignment of the principal axis of p'fp and MMt. 
Notably, the deterministic eddy damping term which has appeared in (^ 

Dni ^ - a< trV^VM >p (lib) 

is a direct generalization of the so-called non-linear eddy damping p^ < tr M^M''' > advanced by Bardina et al. 
p6[ and reduces to it for isotropic tetrads pp'^ = Ip^. In this limit Dni — > — tr{s'^ — O • s • ), where s and fl are 
respectively the symmetric and antisymmetric parts of M. Thus the energy transfer down scale is due to negative 
strain skewness or positive enstrophy production |26j^ (i.e. vortex stretching). We can define the energy flux by 



averaging over the fixed p^ shells. Let R = \J(? and = p'^V°' / R denote the longitudinal velocity, then 

e = -dn (Vn trVV^) ^ + ^ C„e (12) 

is balanced by eddy damping e = Dni + C±Dd- Here (...)^ denotes an additional average over p^ = R^ shell. 

Below we will often think of the diffusive contributions as being small compared to the non-linear interactions on 
current scale: that is, we shall assume C,,, C± <C 1 and treat them as (a singular) perturbation of the deterministic 
dynamics. Another tractable and perhaps physically more plausible limit is C± ^ 1. 



III. DETERMINISTIC DYNAMICS AND THE RESTRICTED EULER MODEL. 

Note that the equation of the form (^^ also governs the Lagrangian evolution of the actual local velocity gradient 
matrix m^h — daVb (we use lower case m to avoid confusion with the coarse grained object) 



dt 



mah 



-dadbP 



(13a) 



as derives from the Euler equation. Leorat 
isotropic contribution to the pressure 



and Vieillefosse pa have considered (13a) retaining only the local and 



dadbP = - ^trm^ 



(13b) 



as a model of vorticity dyna mics and o bserved that (13a, 13b) leads to a finite time singularity with ||m||'^(t, — t)~^. 
The dynamics governed by (13a, 13b ) - the "Restricted Euler dynamics", to use Cantwell's terminology |^5| lies 
entirely in the two-dimensional phase space defined by the two invariants |^,^ triri^ and trxn'^. This reduction 
stems from the 5*^(3) invariance, m —>■ gmg^^ with g being an arbitrary 3x3 matrix, which allows one to bring 
m(t) to diagonal form A{t) by a time independent similarity transformation m(t) = \JA{t)\J~^ . There is yet one 
more independent constant of motion found by Viellefosse [Q: the "discriminant" D = 3{tr (m''))^ — ^{tr (m^))^ = 
— (Al — A2)^(A2 — A3)^(A3 — Ai)^, where Xi{t) are the (in general complex) eigenvalues of m{t). The RE dynamics thus 
reduces to ID flow, i.e. it is integrable! Figure 2 shows the flow in the 2D phase plane of the Cantwell's invariants |3q] 



ktr m^, R 



tr m and the finite time singularity corresponds to the i? — > cx), Q 



Q- 2 

approaching the D = separatrix. 

Along the Z? = separatrix the flow is particularly simple: 

\{t) 
m(i) = I A(<) 

-2A(t) 



-00 asymptotically 



(14) 
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with X(t) = A(0)/(1 — tX{0)) making the finite time singularity at = 1/A(0) exphcit. 

The region above the separatrix, D > 0, is elhptic: the eigenvalues of the velocity gradient eigenvalues become 
complex and the Lagrangian trajectories are rotating; the region below the separatrix, _D < 0, is hyperbolic: the 
eigenvalues are real and the trajectories are strain dominated. These topological aspects of RE dynamics were 
emphasized by Blackburn et al . 

As a model of finite time singularity RE solutions were rejected on the reasonable ground that if considered 
as global solutions of Euler equations these do not satisfy sensible boundary conditions and have unbounded energy. 
Ashurst et al [|j however noted that the statistics of the vorticity/strain alignment observed in the DNS of Navier- 
Stokes may be qualitatively understood in terms of RE. Subsequently Cantwell and coworkers p5|-p7t proceeded to 
investigate the DNS generated statistics of R, Q invariants and observed that the probability distribution function 
(PDF) of R, Q exhibits a pronounced tail along the Viellefosse D = asymptote as can be seen on Fig. 3. These 
two observations suggest that despite the draconian local and isotropic approximation to pressure and the unphysical 
finite time singularity the RE dynamics does capture certain statistical features of the physical flow. 

The deterministic part of the Lagrangian tetrad dynamics defined in Section 2 generalizes RE by reinterpreting the 
velocity gradient tensor as a coarse-grained field defined over the tetrad p and completing the Lagrangian picture by 
adding the dynamical equation for p{t) . The p field introduces the measure of current length scale and the dependence 



on the history of the strain which controls the "shape" of the tetrad. The p-dynamics (3b) is coupled to M via the 
anisotropy tensor 11. For an isotropic tetrad (i.e. regular te trahedro n) Hat = Sab/S and the M-dynamics (the left 
hand side of (pa|)) reduces instantaneously to RE equation (13a, 13b| ). In the next instant however the tetrad will 



become distorted through the action of the volume preserving M following the p-dynamics equation (left hand side 
of (|3b|)) and the trajectory will come out of the RE plane. Its evolution will depart from RE as the anisotropy 
increases and at some point the growth of ||M|| will be cut off. This is most easily seen for the D — Vieillefosse 
line. The dynamics of A (see (|^) becomes A = {6q~^ — 1) A^, with 11 =diag{l, 1, q— 2}/q where q evolves according 
to q = 6{q — 2)A. The isotropic tetrad corresponds to q — 3. Starting from isotropy and A > both X{t) and q{t) 
grow. The growth of q corresponds to the contraction of one of the principle axis of the p^p tensor as the tetrad is 
fiattened in a pancake. However when q > 6 the growth of A reverses. Thus, anisotropy caused by the stretching of 
the tetrad cuts off the Viellefosse finite time singularity [Q. The modified RE dynamics however retains the initial 
growth of M with two expanding and one contracting strain directions and the consequent deformation of the tetrad 
into a pancake or ribbon. This process is the fundamental step of energy transfer |^,^. In the next Section 
we will see that the Viellefosse tail (large i? > 0, Q < region) of the Cantwell PDF on Fig. 3 which is generated 
through this process, indeed corresponds to large negative strain skewness associated with the energy transfer [ p6| . 
Another retained aspect of the RE dynamics is the evolution of the vorticity/strain alignment from configurations 
where vorticity is parallel to a-strain (i.e. the fast stretching direction) to configuration where vorticity is aligned with 
the intermediate, /3-strain as observed numerically ■ The new feature of the modified model is that while in 

the isotropic RE all 2D configurations of M evolve into 3D, (e.g. Mab = ^abz will in the next instant acquire, due 
to low local pressure, a contracting component of strain acting along the z-direction which will act to destroy Wz), the 
new anisotropic model allows the 2D configurations of M to persist provided that the 11 tensor is rank two {p^ p rank 
one) which correspond to filament like tetrads. Note that both intense vorticity and quasi-one dimensional tetrads 
will be produced by the action of strain with one stretching and two contracting directions ir[s'^] > 0, thus there 
potentially is a chance of describing vortex " worms" ||7|,p| JTl]Jl^ . We shall return to the discussion of the kinematics 
of energy transfer and vortex stretching in Section 5. 

How can one compare the deterministic tetrad dynamics model with the real Navier-Stokes dynamics ? The 
relevant empirical object is the average dM/dt and dp/ dt conditioned on M, p but to simplify matters we will restrict 

to isotropic tetrads and examine the flow in the Q^R phase space generated by the conditional averages iR\R,Q 



and \ Q\R, Qj. The latter were obtained by a DNS of Navier-Stokes. 

Briefly, the Navier-Stokes equations are integrated by a standard pseudo-spectral algorithm. Our code is fully 
de-aliased. We used up to (128)^ collocation grid points, and the effective resolution was maintained to be higher 
than fcmax'? > 1-4, where kmax is the highest wave- vector in the simulation, and rjK the Kolmogorov length scale 
{r]K = {v^/ey/'^). Statistics were accumulated for at least 3 eddy turnover times. In the following, we present 
our results for a Taylor scale Re\ = 85. Our investigation of the influence of the Reynolds number in the range 
20 < Re\ < 85 did not reveal any major qualitative change of the statistics presented here. 

Fig. 4a,b,c show the streamlines in the (i?, Q) plane, reconstructed from the conditional averages of (R\R, Q) and 



Q\R, Qj computed numerically for three different p^ . The latter were increasing from the dissipation range to large 

scale. For isotropic p our M-dynamics is instantaneously tangent to the RE and therefore the empirical flows can 
be compared with Fig. 2. Remarkably, while there are signiflcant deviations in the topology of the flow for p in the 
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dissipative range (Fig. 4a), the instantaneous flow for large scale p is surprisingly close to RE. The deviations at small 
scales are presumably due to the viscous effects. The conditional flow for \p\ > lOrjk can be fltted by the modifled 
RE: 

— ^ (a-l)(M^ -UTrM^) (15) 
at 

with a decreasing with increasing \p\/'I]k from .8 to 0. For the reasons related to the energy transfer, discussed in 
Section 2, we believe that a should be constant in the inertial range. The continuous scale dependence observed 
in the flt to DNS however is not unexpected, because the inertial range at the accessible Re is quite limited. On 
the other hand the approximate validity of (^5|) as a description of the coarse-grained Lagrangian evolution is quite 
encouraging. It would be important to extend the comparison of the deterministic dynamics (|l5| ) with the numerical 
simulation for anisotropic tetrads, however in that case the SL(3) invariance of ( p^ ) is lost and in addition to the R, Q 
invariants the time derivative must be conditioned on the vorticity, which makes the computation more demanding 
statistically. It would also be important to investigate systematically the deviations of the conditional flow from (|l^) : 
these are expected to arise from the possible additional deterministic terms in ( p^ ) (e.g. 7M) as well as the stochastic 
dynamics. Much further work is required in this direction. 



IV. LAGRANGIAN PATH INTEGRALS AND THE SEMICLASSICAL APPROXIMATION. 



Let us now explore the statistical properties of the coarse-grained Af-field on the tetrad p. The probability distribu- 
tion P(M, p) is governed by the the Fokker-Planck equation (0,^) but requires specification of an additional boundary 
condition. Since the PDF of velocity is known to be Gaussian on the integral scale we shall fix ||5^ . 



P(M,p)|p2=i2 -exp 



trMMt 



(eX-2)2/3_ 

To impose the integral scale boundary condition one may use a generalization of Green's theorem 
P{M',p')= JdMjdp [F(M,p)Ltgt(M,p|M',p') - LP(M, p)gt(M, p|M', p')] 

ir(dpMV^) P(M,p)gt(M,p|M',p') 



dM 



where denotes the adjoint operator which governs evolution backward in time (obtained by M 
gt(M, p|M', p') = L^~^ its static Green function. 

The static Green's function g^ is computed via the Lagrangian Green's function (0): 



gt(M,p|MV' 



dTGT{M,p\M'p' 



which has an intuitively appealing path integral representation: 



G^t{M,p\M'p') = 



M(a)=M 



DM 



M(-T)=M' 



p(0)=p 



Dp exp[-S{{M,p})] 



P{-T)=p' 



(16) 



(17) 
-M) and 

(18) 
(19) 



summing over all possible paths connecting initial M', p' at time —T with the final M, p at time weighted with the 
action 



S 



dt 



|M - (q - 1) (M2 - mrM2 



C^ep-^ 

tr[{p ~ pM.) (Cl\l - ppt) + C-Vp^) (p - pM)t] 
p^^/trMM.'^ 



(20) 



where = trXX^, p = p/\\p\\ and C|| 0, as assumed before. 
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This path integral form invites a semi-classical approximation 1^,^,^ which estimates the integral via the saddle 
point Gt(M, p|M', p') ~ea;p[— S'c(M, p|M', p')] given by the minimal action Sc along the "classical" trajectories 
connecting the prescribed initial and final points (in time T) and obeying the Euler-Lagrange variational equations. 
Moreover, for each final point (M, p) there exists a unique 5* = trajectory governed by the deterministic part of 
Lagrangian dynamics ( ^ , ph| ) which picks out the Lagrangian preimage M' = M(M, p, — T), p' = /o(M, p, —T) as an 
initial condition. If the small scale generated stochastic component of the dynamics were small Cj_, ^ these 
deterministic Lagrangian trajectories would control the Green's function. Since the probability is constant along 
the zero action trajectory the PDF of the final M, p is determined by the probability of its Lagrangian preimage 
M' = M(M, p, — T) at the integral scale where the PDF is assumed to be Gaussian. Crude as this zero action 
approximation is, it is the natural zeroth order calculation and will provide some physically interesting insights as we 
shall see below. The full semi-classical analysis will be deferred to a forthcoming publication. 



V. PROBABILITY DISTRIBUTION FUNCTIONS AND STATISTICAL GEOMETRY. 



To make contact with the numerical results we shall use the "poor man's" zero action approximation introduced in 
the previous section and according to which the probability of given M observed on a tetrad p in the inertial range 
is equal to the probability of its integral scale preimage. The latter is found by integrating the deterministic part of 
the equations of motion (which generates zero action trajectories) from the observation point backward in time: 

= - (a - 1)(M2 - mrM^) (21a) 



= _gM - M^gt - /3\/trMMt (g - ^ tr g) 



(21b) 



where g = p^ p. The (i term has been added to reintroduce the isotropisation effect due to the transverse small 
scale fluctuations already in the deterministic approximation. This may be thought of as a Mean Field treatment 
of the Cj_ term in (|^) which is physically more appropriate than the formal C^, C„ — > deterministic limit of ^ 



until irg = p reaches the integral scale, yielding for the PDF: 



Dimensionless constants a and (3 will serve as model parameters. Equations ( pla 21bD will be integrated back in time 



P(M, g) ~ exp 



Mi(M, g)M'2(M, g) 



(eL- 



-2^2/3 



(22) 



with M^(M, g) ( and the time of flight T) fixed implicitly by tr [g(M, g, -T)] = L"^. To the extent that the non- 
trivial PDF in this approximation arises as a non-linear mapping of the initially Gaussian variables, our construction 
here is reminiscent of Kraichnan's " Mapping Closure" Pp5[ . 

Let us now construct the PDF of the i?, Q invariants for the isotropic tetrad of radius r. It is convenient to consider 
the elliptic D > Q and the hyperbolic D <{) regions separately and use different parametcrizations of the M matrix: 

a) for D > 



A Ae'' UJ2 

M= I -I\e-^ A Loi 
-2A 



(23a) 



and b) for _D < 



A + A W2 
M = I A - A tJi 
-2A 



(23b) 



where u>a refers to vorticity and in (23a) lo^ 
and h) R= 2X{\^ - A^) and Q = - 



EE 2Acosh7. The invariants are a) i? = 2\{\^ + A^) and Q = A^ - 

3A^ and the strain tensor is s = (M + M^)/2. It is straightforward to 

numerically integrate (21a, 21b ) starting with given M and g — r^l in time until ■^Trg — at which point P(M, g) 
is assigned via (p2). However, to determine Pr{R,Q) one must integrate P(M) over a) 7, tJi.2 and b) over a;. (Note 
the Jacobian: J dMS{trM) = J dRdQduj = J d\dAduj^/\D\ and dAduJi = 2A7 sinh [7] d7dA). 
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The task is simplified within the saddle approximation where the integration is reduced to minimization K] of 



logP(M.^ r) with respect to the integration variables which we carry out numerically via an "amoeba" algorithm 



57 



Over the whole R, Q plane we find that the saddle point is at a) 7 = a;i^2 = and h) uj — 0. In addition for the special 



case of D = 0, when A is the only non-zero parameter in (23a, 23t) the trajectory and the P(M, g) can be computed 



analytically (see Appendix A). The resulting distribution (for different r) is displayed on Fig. 5a, b. P{Q,R) exhibits 
a long (but Gaussian ^exp{~aX^r'^''+)) ridge - the " Vieillefosse tail" - along the D = line in the i? > 0, Q < 
quadrant and a valley of low probability approaching the origin from the R > side. This structure appears because 
the backward in time trajectory of the point on the D = separatrix converges to the origin and maps to a highly 
probable, small M' integral scale preimage, whereas the trajectories originating in the low probability gulf in their 
time reversed dynamics are swept westward past the origin into the improbably large M' region. The PDF Pr{R, Q) 
evolves continuously with decreasing r/L away from the Gaussian (which appears non-trivial in the i?, Q variables) 
shown for comparison on Fig. 6. The appearance of the PDF tail along D = 0, R > Q and the the trend of r/L 
dependence are reminiscent of those for the PDF obtained from the Navier-Stokes. Yet, both high probability ridge 
and the low probability valley found in the present deterministic approximation are strongly exaggerated. There is 
a good reason to expect that the effect of the fluctuations will be strong in these regions: the narrow ridge should 
be largely "washed out" as the asymptotic behavior of the tails and low probability regions is clearly dominated by 
fluctuations. Yet again, the rather complex i?, Q dependence of the PDF and the crude qualitative similarity of the 
model PDF with the results of the DNS merits a detailed discussion of the underlying kinematics and dynamics. 

Let us compute the distribution in i?, Q plane of the average enstrophy uj"^ and enstrophy production uj ■ s ■ uj — 
iTr(]V[ -t- M'I')(M — M'l')^ which measures the rate of vortex stretching ||2^. The enstrophy density is defined by 
e(i?, Q) = J dLuuj^Pr{R, Q,oj) but in order to save computer time will only be evaluated in the saddle approximation 
by varying the integrand w.r.t. uj. For all of the hyperbolic region the saddle of the integral occurs at a non-zero 
values of 7, uj parameters (21^,21b). The result is presented in Fig. 7. We see that the average enstrophy peaks at 



i? = and small positive Q. This is explained by noting that the conditional average Q) grows like Q (at 

least for i? « and Q > 0) because Q = {\uj^ - Tr s^)/2, while P(i?, Q) falls off. The average enstrophy production 
a{R,Q) = J dujuj ■ s ■ uj Pr{R,Q,u}) is also dominated by a non-trivial saddle and has the R,Q dependence shown 
in Fig. 8. We observe that enstrophy is produced predominantly in the upper left quadrant of the R, Q plane, it is 
(partially) destroyed in the upper right quadrant and there is weak vortex stretching in the D = tail. The tail region 
is not well resolved on Fig. 8 and the positive enstrophy production is confined to the narrow strip delineated by the 
zero contours. This strip becomes progressively more and more narrow as one looks at smaller scales, r, (which is 
why we limited our figures to rather large r/L,) however, the existence of positive vortex stretching in the tail domain 
was verified analytically via the arguments presented later in this section (and in Appendix A). Next we compute the 
average strain skewness, S3 = J duj tr s"^ Pr{R, Q) , which is the object associated with energy transfer |2^, and find 
that it is strongly localized in the D — Vieillefosse tail as shown in Fig. 9. (Curiously, we find that in the elliptic 
region, Q > 0, the strain skewness changes sign three times - a fact that can be understood on the basis of the structure 
of the M tensor given by (25).) Note that the energy transfer term in (|^) is actually Tr M^M^ = Trs^ — juj • s • w so 
that vortex stretching also contributes to the energy fiux. Q Vortex stretching however does not dominate the energy 
transfer: while most of the enstrophy production is in the upper left quadrangle. Fig. 8, the energy flux distribution 
F{R,Q) EE -a J dujTrM^M'<Pr{R,Q,uj) is localized in the D = tail where the negative Trs^ lives, see Fig. 10. 
Thus we arrive at the conclusion that from the energy transfer point of view, the active regions are dominated by the 
strain and not vorticity [ p^p5[ . Also of interest is the appearance of distinct regions of weak negative energy transfer: 
e.g. the R,Q > quadrangle where both a and the energy flux are negative [Q. 

Remarkably, as seen from the comparison of Fig. 7-10 and Figs. 11-14 the DNS exhibits a rather similar distribution 
of average enstrophy, enstrophy production and strain skewness in the Q, R plane. Both DNS and the model have 
positive enstrophy production in the upper left quadrant and in the D = 0, R > tail, and negative cr in the 
upper right quadrangle; both have strain skewness strongly confined to the D = 0, R > tail, and both exhibit the 
dominance of the strain skewness in the energy transfer. Furthermore, there is a clear correspondence of the positive 
and negative regions of <rs^ (Fig. 9 and 13) and energy flux (Fig. 10 and 14). 

Much of this behavior can be understood by considering the "statistical geometry" of the flow (see Ashurst et al 
Constantin ||l^] and Tsinober et al |2^] for excellent discussions) and is inherited from the RE dynamics. Most of 
the vortex stretching occurs for R < 0, Q > where the vorticity is aligned with the large positive eigenvector of the 
strain; this is immediately evident from the structure of the M matrix at the saddle point controlling cr(R, Q): 




= -A A (24) 

\ -2A y 

where p^A < 0, since the off-diagonal elements of (24) define the vorticity along the direction with the strain 
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eigenvalue — 2A > 0. The tetrad dynamics then takes the M field into the i? « 0, Q > region where A and hence the 
stretching rate a vanish before changing sign. The magnitude of vorticity reaches a maximum. In the peak vorticity 
region M tensor has the form: 



Me = -Ae-^ A (25) 





with A ^ A and 7 ^ 1, so that the flow is essentially two dimensional with vorticity aligned with the nearly neutral 
strain direction. 

The alignment of intense vorticity with the intermediate strain axis is well known [^j^j^,^,^^ and its reappearance 
in the model is encouraging. The fate of the typical ID vorticity is to implode under the action of the contracting 
strain (brought about by the low pressure of the vortex) pof . The energy is transfered to larger scale. 

To conclude we examine the geometry of M in the _D = tail. Here the predominant M configuration is 

A w 

= I A I (26) 



-2A 



which leads to the strain eigenvalues 



=A; -^V'+9A^-2 



and s-y = —Sa — sp. The vorticity is aligned with the intermediate strain axis, which is stretching since the sp > 
consistent with the negative strain skewness Trs^ ~ SsaSps^ The normalized /^-strain {^1'^) ^ function 

of A along the D = line increases as ^^|A^ goes down with increasing A as it is shown below in the Section and 

confirmed by Fig.l5a,b for the DNS. To the extent that Pr{R,Q) is peaked at i? = Q = and hence at A = the 
most probable configurations have si ^ S2 ^ and are close to plane shear or 2D vorticity. Note that these most 
probable configurations have Tr « or Trs^ « which means approximate local homogeneity - a much 
stronger statement than the homogeneity on average, i.e. (TrM2) =0! This must have serious implications for the 
pressure. 

The exagerated non-Gaussianity (i.e. Vieillefosse tail) of the deterministic model aside, there are also other quali- 
tative differences between the results of the model and the DNS: e.g. the vorticity density in the DNS peaks at the 
origin and is skewed in the upper half plane, towards R < 0. This effect is even stronger for the vorticity distribution 
measured in the dissipative range: an asymmetry also visible in P{R, Q) shown on Fig. 3a. This can be understood 
as a result of enstrophy dissipation by viscosity which shifts the locus of null enstrophy growth and hence maximal 
enstrophy to configurations where of vorticity is still aligned with positive strain. 

Let us further explore the results of the deterministic approximation. The exact solution on the D = Q line (see 
Appendix A) exhibits strong asymmetry: even though the PDF is Gaussian along the D — line both for R > and 
R < 0, lnP{X) ~ — A^(r/L)^(^~"^''±('') the characteristic exponents on the two sides are different, h+{[3) > h-{[3). 
It is the /i_ exponent which controls the behavior of low moments including the 3d and in order to impose constant 
energy flux we must require (1 — a)h^{(3) = 2/3 which forces K41 scaling on the "head" of the PDF in Fig. 5a, b. The 
constant flux condition thus fixes a specific relation (determined in the Appendix A) between the model parameters: 
as /3 increases from to « .3, a decreases from 1/3 to 0. 

However, while the scaling in the bulk of the PDF is K41, the scaling in the Vieillefosse tail is anomalous. Introducing 
rescaled variable A = Ara, A = Ars, and uj — Lor^ we find that the A dependence of the action in the hyperbolic 
vicinity of the tail (which determines its width in the R, Q plane) has the asymptotic form: 



S'+(A, A,w = 0) - AV^Tt + AV-^Va --r-''-T+ (28) 




where the scaling function /a(0) = const and fA{x)^x'' for x S> 1 with 77 = 7+C/(2 — C). Similarly, the vorticity 
distribution on the D = {R > 0) line is governed asymptotically by 

5+(A, A = 0,Cj) - A2r2T+ + Cj^'^'^' J^^^-V-T+j (29) 
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where fuj{x) ~ for a; <C 1 and /i^(cc) ^ ^ for 2; ^ 1 with rj' = 7+ (3 — S)/{5 — 1). The anomalous exponents 
depend on /?: 7+ = 2(1 - a)(/i+(/3) - h-{P)) varies from at /3 = 0, (a = 1/3) to 1/2 at /3 « .3, a = 0. Over the 
same range 77 varies from 1 to 3/2. The exponents as a function of /3 are tabulated in the Appendix A 
It follows that the conditional enstrophy in the Vieillefosse tail behaves for large A as: 

< u^\~X (A)^ r^+2V (30) 

implying an interesting non-trivial scaling relation between the strain and vorticity of the velocity gradient sheets 
associated with the high energy transfer regions. 

We can also compute the contribution of the D 0, i? > tail to the monts of the velocity gradient. E.g. we 
compute 

(a")^ ^ = N JdX^ J dA^ J dc:;A"e-^+(^'^'+'^''"''^^'-""') ^ r2^+3V-(«+3)7+ (31) 

where the normalization factor ~ 1 because it is dominated by the "head" of the PDF with K41 scaling, already 
factorized explicitly. Thus we conclude that the tail contribution is important only for sufhciently large n, when 
n > ric = (2?7 + 7^ — 3 and the anomalous scaling becomes dominant over the normal contribution of the 

bulk of the PDF. We find that in the limit of /3 and a ^ 1/3, nc —> 00 and the shape of the PDF, although 
non-Gaussian, becomes independent of the scale. This limit recovers the Kolmogorov 41 theory. The intermittency 
effect is maximized as /3 — > .3 and a ^ 0, where anomalous scaling appears for n > 6. 

It is clear that although the PDF found in the present deterministic approximation of the model is on the whole 
far from Gaussian (imposed on integral scale) the approximation underestimates the intermittency effects. E.g. it 
predicts no deviations from K41 in the vorticity dominated Q > region, in contrast to the DNS result which in that 
region exhibits r-dependence associated with the development of an exponential (or sub-exponential) tail on smaller 
scales. It is equally clear that the description of such a tail is beyond the currently employed approximation. After all, 
the asymptotic behavior of the PDF, i.e. the statistics of large fluctuations, is dominated by fluctuations. One does 
however expect to find the exponential asymptotics ]4l| ] once the fluctuations are accounted for via a proper semi- 
classical calculation. There thus appears to be two mechanisms contributing to the intermittency: 1) the deterministic 
non-linear interactions in the energy transfer region as indicated by our present calculation and 2) the effect of small 
scale fluctuations responsible for the exponential asymptotics of the PDF. The latter mechanism is analogous to the 
one responsible for the intermittency of the Passive Scalar |l^,^,^,|8| . 



VI. CONCLUSIONS. 



In the preceding Sections we have introduced and began to analyze the phenomenological model of inertial scale 
velocity fluctuations defined through the velocity gradient tensor coarse-grained over a region specified by a tetrad 
of points. The dynamics of this field was decomposed into the non-linear deterministic component representing local 
same-scale interactions and a Gaussian stochastic component with embedded K41 scaling representing interactions 
non-local in space and the incoherent contribution of the velocity fluctuations from scales smaller than that of the 
tetrad. The deterministic component is closely related to the RE dynamics of Vieillefosse model. The latter, 
although marred by the unphysical finite time singularity, ]6l|]has been an appealing candidate description for the 
dynamics of the local velocity gradient [p|, p^|j39| ] . A novel aspect of our model is the elevation of the velocity gradient 
dynamics to the coarse-grained level and the introduction of the tetrad tensor (p) dynamics which explicitly introduces 
the scale and the measure of anisotropy generated through the strain induced distortion of the Lagrangian volume. 
This allows to construct an anisotropic model of the coarse grained pressure Hessian which eliminates the finite time 
singularity from the deterministic dynamics by suppressing the work done by the pressure on the distorted fiuid 
element. Furthermore, the explicit appearance of the current scale allows to build in K41 spectrum in the stochastic 
component p7[ of the dynamics. An important reason for working with the coarse grained field is that in the dynamics 
of the purely inertial range fields the direct contribution of the viscosity can be neglected |6^. Instead, the ultimately 
viscous dissipation is incorporated through the effect of the incoherent small scales acting through the eddy damping 
D„i + D± terms in the model. The deterministic dynamics plays the key role in transferring energy down scale. This 
transfer occurs due to the volume preserving distortion of the fluid element which leads to the reduction of at least 
one of its principle axis. 

Our heuristic derivation of the model was fortified by the numerical test of ( ^ ) through the construction of the 
conditional fiow for the invariants of M. The deterministic part of (3a) appears to be quite close to the conditional 
flow at least for isotropic tetrads p = Ir. It will be important to extend the numerical study to anisotropic p and to 
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find a way of examining the stochastic contribution to the dynamics. The vahdity of the local approximation of the 
pressure Hessian jl^ and the neglect of non-local correlations remain the key issues. 

The tetrad model respects K41 scaling (M — > foM, p b^^/'^p), both on the operator (0) and the integral scale 
boundary condition (|l6| ) levels. Yet, as illustrated by our crude deterministic approximation the resulting PDF has 
non-trivial, model parameter dependent anomalous scaling. The Kolmogorov scaling has to be re-imposed on the 
level of the 3d moment by choosing the parameters so as to make the energy flux scale independent. Except for 
one particular limit (/3 — + 0, a — + 1/3) where the PDF has K41 behavior, high moments exhibit anomalous scaling. 
Remarkably, while this anomalous scaling has little to do with the Kolmogorov-Obukhov arguments |26|Jll| it does 
originate in the domain of energy transfer: the Vieillefosse tail. 

The deterministic approximation employed in Section 5 grossly overestimates the extent of the Vieillefosse tail in 
P((5, R) but it does bear resemblance to the PDF observed in the DNS. It also generates plausible distributions 
of enstrophy, enstrophy production and strain skewness. The present analysis provides a clear dissection of the 
high enstrophy and the high enstrophy production regions, identifies the difference in the vorticity- strain alignment 
in the two regions and exhibits their dynamical connection. As explained in Section 3, much of this sensible 
phenomenology was inherited from the RE dynamics The elimination of the Vieillefosse finite time singularity 
was however essential in order to have a model with stationary statistics. 

The energy transfer for isotropic p configurations occurs via non- hnear eddy damping @,|§ irM^Mt term (|), 
which combines contributions of strain skewness and vortex stretching. Within our present crude deterministic 
approximation, the energy transfer occurs largely in the Viellefosse tail and is due to large negative trs^ . In this 
region strain has two positive eigenvalues and Lagrangian volumes are deformed into pancakes or ribbons, i.e. this 
is the region of sheet formation In contrast, the vortex filaments are generated in the Q > 0, i? < quadrant 
where enstrophy production is peaked. This region does not contribute as much to the energy fiux. Furthermore, 
the maximal vorticity region is characterized by nearly 2D configurations and does not contribute at all - a notion 
consistent with recent numerical results p^ . Yet, whereas we are optimistic about correct description of the high 
—trs^ tail in the model, the correct description of the high a region may be more difficult because of the importance 
of long-range strain in vortex stretching. 

In order to enable a more quantitative comparison of the model and the DNS the treatment of the model must 
include the effect of the fluctuations which control the asymptotic behavior of the PDF. Our present, deterministic, 
approximation overestimates the PDF in the narrow Vieillefosse tail, while underestimating the asymptotic behavior 
of the PDF elsewhere in the Q, R plane. Both effects are due to the neglect of the fluctuations. Since the relevant 
Fokker-Planck equation lives in the (P — 1 -\- d{d + l)/2 = 14 dimensional M, p space the direct numerical approach 
seems out of question. However one may hope to make progress with the semi-classical analysis of the path integral 
representation (|2^) along the lines of |^,^ . This approach is valid for the calculation of the PDF tails and associated 
anomalous scaling. We expect that the fluctuation effects will change the Gaussian decay of the PDF tails found in 
the deterministic approximation to the exponential, e.g.: A~" e"'^^'' . 

An interesting simplification appears in the limit of strong transverse diffusion Cj_ > > 1 which corresponds to the 
physically sensible regime of strong re-isotropisation of the tetrads arising from the action of incoherent small scale 
fluctuations. In that limit the PDF becomes nearly uniform over the p^ — const shells and can be projected onto the 
s and d representations of 5*0(9) acting on pf: 

P{M,p) « *(M,p2) + - ^P'S'^'S,, ) $^''(M,p2). (32) 



Curiously, this nearly isotropic approximation doubles as a d >> 1 expansion [p8| , [L9[|29| ] because the d-wave mode 
of SO{d^) is suppressed by the l{l -I- c?^ — 2) total angular momentum factor with I = 2. This expansion allows a 
considerable simpliflcation of (|^) and (pO|). 

Several related models are worth mentioning. The present model appears to be the real space, Lagrangian coun- 
terpart of the momentum space, Eulerian "shell" model proposed by Siggia ||7^. Yet, the geometrical and statistical 
implications of that model have not been fully explored and we do not at present understand the relative merits of 
the assumptions involved in the two models. The logic which led to (^ 3b) can and has been applied to the Passive 
scalar problem, in which case one would replace M by a Gaussian random field. This would lead to a PS model of 
the type considered in p0|j 7l ] closely related to the Kraichnan's ||2^ model. It would be interesting also to explore 
whether the 2D version of (|3a,3b) could provide a sensible description of 2D turbulence where the physics is very 
different. In that case the deterministic left hand side of ( pa] ) would vanish for isotropic "triads" . We do not presently 
know if the model can generate the inverse cascade. Finally, if the tetrad model provides a reasonable description of 
turbulent fiuctuations in the homogeneous case, it will be interesting to attempt to generalize it to the inhomogeneous 
and anisotropic case: e.g. one could study the statistics of tetrads moving away from the wall in the boundary layer. 

Last but not least, two remarks concerning the relation with the experiment. The traditional approach to turbulence 
and most of the existing data involves velocity difference at two points. This two point statistics can be extracted 
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from the tetrad statistics by averaging over "unobserved" variables P(i;, pi) = / dp2 J dp^ J dM5(p"Mab — W6)P(M, p). 
Conversely, it would be very interesting to study tetrad statistics experimentally. For that we need a velocity mea- 
surement at 4 points in the inertial range, which can be obtained from 3 (crossed wire) probes (say a fixed probe at 
the origin and two movable probes at {x,y,z) = (0, y/3/2r,±l/2r)in a flow with < Vx >^ 0. The time lag (on the 
static probe) then can be used to provide the 4th measurement jT^J?^. The PDF coarse grained velocity gradient 
may then be obtained. Quite independently of the predictions of the present model, the (i?, Q)-plane density of the 
2nd and 3d order invariants (e.g. Figs. 8-14) are quite illuminating in dissecting the role of vorticity and strain in 
intermittency and the possible difference in their anomalous scaling. 

We conclude finally that despite its relative simplicity the tetrad model is surprisingly rich in physics, offering an 
insight both into the geometry and dynamics as well as the statistical and scaling properties of the inertial range fields. 
Many non-trivial statistical objects can be calculated in terms of the 3 parameters of the present model. Further 
work and detailed comparison with the DNS and the experiment will establish the degree of success and failure of 
this model. 
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APPENDIX A: SCALING IN THE DETERMINISTIC LIMIT 



Consider Lagrangian dynamics of M and g described by (21a, 21b ) for the following class of matrices forming a 
subspace in the Elliptic region {D > 0) 



M = 

We arrive at the system of equations 

din [A] 

Xdt 
d In [/i] 
Xdt 
dz 
Xdt 
din [G] 
Xdt 



X 


A 





-A 


A 











-2A 



X 
\ X 
y 



(l-a) 
{1-a) 



/i z 



z + 2 
z + 8 + p'^ 
z + 2 



-6z + s.gn[A]|V6 + V(l + 2z)(l - z), 

4(1-^) 
l + 2z ' 



(Al) 

(A2) 
(A3) 
(A4) 
(A5) 



with /i = A/A, z = x/y (z(0) = 1), G = tr[g] ~ 2x + y. This system is integrable in quadratures if either f3 — or 
A(0) — (and therefore, A{t) = 0, for any other t). We shall now calculate the effective action 5 = — In [P(M, g)](see 
(pa)) for these two cases. 



1. Elliptic region, (3 = 0. 



Integration of the system ( [A2| - |A5| ) gives 

A'2: 



2r6 



X X'^-A^3°'-^ 



(A6) 



r^Y " /l-a 4-a -2r^ 



2F1 



l—a 4— a 



-,a. 



3 ' 



A' = A 



2 \ 



(A7) 



13 



where the "' notation stands to mark the (— T) preimage. These expressions allow to rewrite the effective action in 
terms of r, A and A, providing T is fixed by tr{p{—T)p{—T)^] = 2x' + r^/a;'^ = 3L^ = 3. There are two solutions for 
x' realized separately depending on the sign of A'. 

First, consider the region with A being positive during the all Lagrangian evolution. At r <C 1, the respective value 
of the effective action is 



A' > 0, 5 ^ 3A^ 



2 \ 



V3 

r 



2(l-a) 



(A8) 



We find that the action becomes infinite at r — *■ in the domain. 

If A(<) is negative (at least at the final stage of the backward in time evolution) and r ^ 1 the effective action is 



A' < 0, 



S 



y2 ^ ^230-1 _^p^ 



1 — a 4 — a 



3r2 
~2~ 



l-a 



(A9) 



The transition region between ( |A8[ ) and (A9) shrinks with r ^ 0. The crossover occurs at the intersection of the 
A = /xA line in the A, A ( Q, R) plane, /i, depends on both a and r and has the following asymptotic form: 



3Q-I |"2fa-11/3p 



/[-3r(a-l)] -2F1 



' rr — 

3 ' ' 3 



3 / V 3 

3i/2-"(i„Q,)r4a-i/[2"+i(l - 4a)(2a + 1)], a< 1/4 



2)] , a> 1/4; 



(AlO) 



Therefore, a sector in the right part oi Q,R plane bounded by the D = Q line from below and the A = fiA one from 
above forms a low probability gulf of the PDF. 



2. Zero discriminant line: D = Q. 



For the 6 — D — line, the integration of ( |A2| - |A5| ) yields 



In 



G' 



In 



3r2 

A_ 

A' 



= 4 



(1 — z)dz 



6i - /3 (1 + 2z) (1 - ys/Ssgn [A] (1 + 25) 



= (!-«) 



(4 - z) dz 



( (i{l + 2i){l-z).j2/isgn[\]^{2 + i) 



(All) 



The dynamics (on D — line) does not change the sign of A so that sgn[X\ is constant. 

For positive A, \{t) is monotonically decreasing while G{t) grows and z{t) decreases. For small enough initial 
G(= 3r2 <C 1), z approaches 2;+, 



f3 - 3\/6 + ^54 - 6/3V6 + 9/32 
4^3 • 



For z' close to one finds 



where 



A' 

— > const * [z — Z- 

A 



c+ = V372 



,C+(l-a) 



G' 



const * 
[z' z+f^ 



4-Z4 



C+ = V24- 



(2 + z+)V54-6/3V6 + 9/32' 

1-Z+ 



(1 + 2z+)v/54-6/3%/6 + 9/32 
We get the following asymptotic behavior for the effective action for Z? = and i?, A > 0: 



(A12) 

(A13) 

(A14) 
(A15) 
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(A16) 



with h+ = C+/C2+. 

If A is negative, X(t) increases in absolute value while z{t) decreases at the initial stage of backward in time evolution. 
The dynamics changes ( r is supposed to be very small) once z crosses 4 and X{t) starts moving towards the origin. 
z keeps growing to approach z_ , 



3^/6 + f3+ V54 + 6V^^ + 9/32 
4^ 



For z close to z_ one finds 



where 



const * \z — z_ 



C2" = V2A 



|C-(l-a) 



-4 



const * r 



(2 + z_)V54 + 6V6/3 + 9/32' 
z_ - 1 



^ (1 + 2z_) v/54 + 6%/6/3 + 9/32 
Finally, we get the following asymptotic behavior for the effective action for Z) = and i?, A < 0: 



(A17) 

(A18) 

(A19) 
(A20) 

(A21) 



with h- = Ci /C2 ■ Note, that < h+. 

Remarkably, the numerical study of the action shows that this scaling found analytically for the D — 0, R < line 
holds everywhere in the Q, R plane except for the Vieillef osse tail, D — 0, R > 0. Hence the main body of the PDF, 
which determines the low moments, scales according to (A21). In order to impose constant energy flux we fix the 
scaling of the low moments (including the 3d) to the K41 value, which requires 



a = 1 - 



1 

3hZ' 



(A22) 



thus relating a and (3 parameters of the model, a decreases with /3 increase from 1/3 at /3 = (an exceptional case 
when tail and the body of the PDF scale the same way) to at /3 « 0.3. We keep /3 as free parameter in the interval 
[0, « 0.3] and calculate the respective values of other exponents numerically Table 1. 
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FIG. 1. Four Lagrangian points forming a "tetrad". The velocities at the four points define the coarse grained velocity 
gradient field. 
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FIG. 3. The PDF of Q*,R* invariants normalized to the variance of strain, = Q/ < > and = R/ < >^/^ ("star" 
denotes normalization) , obtained from DNS at 7?^ ~ 85 measured at different length scales: a) dissipation range p = 2ri h) low 
end of the inertial range p = 8r;, and c) upper end of the inertial range p — L/2. The isoprobability contours are logarithmically 
spaced, and are separated by factors of 10. The dashed line corresponds to zero discriminant. 
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FIG. 5. PDF of Q,, JZ* invariants (normalized as in Fig.3) calculated for the tetrad model in the deterministic approximation 
a) p/L = .2; b) p/L = .5. 
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FIG. 7. Enstrophy density in R, Q plane for the tetrad model at p/L = .5 
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FIG. 8. Enstrophy production density in R,Q plane for the tetrad model at p/L = .5 
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FIG. 9. Strain skewness density in R, Q plane for the tetrad model at p/L = .5 
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FIG. 10. Energy flux density in R, Q plane for the tetrad model at p/L = .5 
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FIG. 11. Enstrophy density in R,Q plane from the DNS, Rx = 85, at p/L = 0.125 (same value as in Fig.3c,4c). The 
enstrophy is normalized by —2 < tr(r2^) >. 
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FIG. 12. Enstrophy production density in R, Q plane from the DNS, R\ = 85, at p/L = 0.125. The enstrophy production is 
normahzed by | < tr{M^ A4^) > |. Solid lines correspond to positive values, dashed lines to negative values. 
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FIG. 14. Energy flux density, tr{M'^M'') in R, Q plane from the DNS, Rx = 85, at p/L = 0.125. The energy flux density is 
normalized by | < tr{M^ M^) > |. The same convention as in Fig. 12 is used. 
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FIG. 15. a) Normalized intermediate eigenvalue of the strain , s« = [6s/3]^''^/|s|, as a function of A along the Vieillefosse 
tail as observed in the DNS, R\ = 85, p = 2rj\ b) same for normalized enstrophy, e* = (a;^|A)/A^, as a function of A. 
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0.25 


0.33 


0.43 


V 


1.11 


1.14 


1.21 


1.28 


1.35 


1.44 


C 


1.95 


1.86 


1.71 


1.66 


1.63 


1.55 


s 


1.33 


1.69 


1.69 


1.85 


1.89 


1.90 


v' 


0.27 


0.22 


0.34 


0.34 


0.42 


0.52 


Tic 


62.98 


28.87 


18.82 


12.88 


9.92 


7.52 



Table 1 Table of exponents for the D = 0, R> tail in the deterministic limit (Appendix A). 
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